ens <- read.delim("mart_export.txt")
u <- duplicated(ens[,1])
ens <- ens[!u,]
p <- ens[,5] == "protein_coding"

x <- dimnames(exp_matrix)[[1]] %in% ens[p,1]
exp_matrix <- log2(exp_matrix[x,]+1)
x <- apply(exp_matrix,1,sum)
exp_matrix <- exp_matrix[x>0,]



cv <- apply(exp_matrix,1,var)/apply(exp_matrix,1,mean)
pickSoftThreshold(t(exp_matrix[cv>0.089,]))
adj_mat <- cor(t(exp_matrix[cv>0.089,]))
a <- adj_mat < 0
adj_mat[a] <- 0
tom_mat <- TOMdist(a^6)
clust <- hclust(as.dist(tom_mat),method="av")
mods <- cutreeDynamic(clust,distM=tom_mat)
ME <- moduleEigengenes(t(exp_matrix[cv>0.089,]),mods)
kME <- cor(t(exp_matrix),ME$eigengenes)

##For plotting MEs
wells <- read.csv("well_names.csv")
a <- strsplit(wells[,1],"_",fixed=TRUE)
b_nc <- rep("",96)
b_geno <- rep("",96)
b_tx <- rep("",96)
b_duration <- rep("",96)
for(i in c(1:96)){
	b_nc[i] <- a[[i]][2]
	b_geno[i] <- a[[i]][3]
	b_tx[i] <- a[[i]][5]
	b_duration[i] <- a[[i]][4]
}
exc <- b_nc == "NC230614" | b_nc == "NC230829"
b_geno_num <- rep(0,96)
b_geno_num[b_geno == "HET"] <- 1
b_geno_num[b_geno == "NULL"] <- 2
b_tx_num <- rep(0,96)
b_tx_num[b_tx == "RAPA"] <- 1
b_dur_num <- rep(0,96)
b_dur_num[b_duration == "Acute"] <- 1
o <- order(b_geno_num[!exc],b_tx_num[!exc],b_dur_num[!exc])
col_geno <- rep("blue",96)
col_geno[b_geno=="HET"] <- "purple"
col_geno[b_geno=="NULL"] <- "red"
d <- rep(0,96)
d[b_duration=="Chronic"] <- 10
barplot(ME[[1]][!exc,"ME2"][o],col=col_geno[!exc][o])
barplot(ME[[1]][!exc,"ME2"][o],density=d[!exc][o],add=TRUE)

##For heatmap
col_geno <- rep("blue",96)
col_geno[b_geno=="HET"] <- "purple"
col_geno[b_geno=="NULL"] <- "red"
col_tx <- rep("cyan",96)
col_tx[b_tx == "RAPA" & b_duration == "Chronic"] <- "green"
col_tx[b_tx == "RAPA" & b_duration == "Acute"] <- "magenta"

a <- cor(exp_matrix[cv>0.089,])
h_cols <- colorRampPalette(c("purple","red","yellow"))(64)

##Repeat WGCNA without extra samples
x <- dimnames(exp_matrix)[[1]] %in% ens[p,1]
exp_matrix <- log2(exp_matrix[x,]+1)
x <- apply(exp_matrix[,!exc],1,sum)
exp_matrix <- exp_matrix[x>1,!exc]

mn <- function(x){
	out <- rep(0,12)
	count <- 1
	for(i in c(0:2)){
		for(j in c(0:1)){
			for(k in c(0:1)){
				out[count] <- mean(x[b_geno_num[!exc] == i & b_tx_num[!exc] == j & b_dur_num[!exc] == k])
				count <- count + 1
			}
		}
	}
	return(out)
}
condition_means <- t(apply(exp_matrix,1,mn))
cv <- apply(condition_means,1,var)/apply(condition_means,1,mean)

pickSoftThreshold(t(exp_matrix[cv>0.018,]))
adj_mat <- cor(t(exp_matrix[cv>0.018,]))
a <- adj_mat < 0
adj_mat[a] <- 0
tom_mat <- TOMdist(a^6)
clust <- hclust(as.dist(tom_mat),method="av")
mods <- cutreeDynamic(clust,distM=tom_mat)
ME <- moduleEigengenes(t(exp_matrix[cv>0.018,]),mods)
kME <- cor(t(exp_matrix),ME$eigengenes)

z <- rep(0,88)
z[b_geno_num[!exc] == 2 & b_tx_num[!exc] == 0] <- 1
z[b_geno_num[!exc] == 2 & b_tx_num[!exc] == 1 & b_dur_num[!exc] == 1] <- 1
mod_correlation <- cor(ME[[1]],z)

make_all_pcs <- function(ME,f){
	n <- dimnames(ME[[1]])[[2]]
	pdf(f)
	for(i in c(1:dim(ME[[1]])[[2]])){
		barplot(ME[[1]][o,i],col=col_geno[!exc][o],main=n[i])
		barplot(ME[[1]][o,i],density=d[!exc][o],add=TRUE)
	}
	dev.off()
}

net <- cbind(dimnames(kME)[[1]],ens[dimnames(kME)[[1]],6],m_names[m],kME[cbind(1:dim(kME)[[1]],m)])

##Mito
mito_overlap <- matrix(0,41,2)
m_names <- unique(network[,3])
for(i in c(1:41)){
	z1 <- network[,3] == m_names[i]
	z2 <- network[z1,2] %in% mito[,3]
	mito_overlap[i,1] <- sum(z2)
	mito_overlap[i,2] <- phyper(sum(z2),sum(z1),18507-sum(z1),1157,lower.tail=FALSE)
}

##FMRP
fmrp_overlap <- matrix(0,41,2)
m_names <- unique(network[,3])
for(i in c(1:41)){
	z1 <- network[,3] == m_names[i]
	z2 <- tolower(as.character(network[z1,2])) %in% tolower(as.character(fmr.d[,2]))
	fmrp_overlap[i,1] <- sum(z2)
	fmrp_overlap[i,2] <- phyper(sum(z2),sum(z1),18507-sum(z1),761,lower.tail=FALSE)
}

##EGR1
x <- egr1[,6] > 10 & egr1[,11] > 0 & egr1[,11] < 5000 & tolower(as.character(egr1[,10])) %in% tolower(as.character(network[,2]))
egr1_targets <- unique(egr1[x,10])

egr1_overlap <- matrix(0,56,2)
m_names <- unique(network[,3])
for(i in c(1:41)){
	z1 <- network[,3] == m_names[i]
	z3 <- network[z1,4] > 0
	z2 <- tolower(as.character(network[z1,2][z3])) %in% tolower(as.character(egr1_targets))
	egr1_overlap[i,1] <- sum(z2)
	egr1_overlap[i,2] <- phyper(sum(z2),sum(z3),17469-sum(z3),length(egr1_targets),lower.tail=FALSE)
}
dimnames(egr1_overlap)[[1]] <- m_names

##Repeat WGCNA without extra samples and outliers
exc <- b_nc == "NC230614" | b_nc == "NC230829" | colnames(exp_matrix) %in% c("D7","H8","F5","A10")
x <- dimnames(exp_matrix)[[1]] %in% ens[p,1]
exp_matrix <- log2(exp_matrix[x,]+1)
x <- apply(exp_matrix[,!exc],1,sum)
exp_matrix <- exp_matrix[x>1,!exc]

cv <- apply(exp_matrix,1,var)/apply(exp_matrix,1,mean)

pickSoftThreshold(t(exp_matrix))
adj_mat <- cor(t(exp_matrix))
a <- adj_mat < 0
adj_mat[a] <- 0
tom_mat <- TOMdist(adj_mat^6)
clust <- hclust(as.dist(tom_mat),method="av")
mods <- cutreeDynamic(clust,distM=tom_mat)
ME <- moduleEigengenes(t(exp_matrix),mods)
kME <- cor(t(exp_matrix),ME$eigengenes)

z <- rep(0,84)
z[b_geno_num[!exc] == 2 & b_tx_num[!exc] == 0] <- 1
z[b_geno_num[!exc] == 2 & b_tx_num[!exc] == 1 & b_dur_num[!exc] == 1] <- 1
mod_correlation <- cor(ME[[1]],z)

make_all_pcs <- function(ME,f){
	n <- dimnames(ME[[1]])[[2]]
	pdf(f)
	for(i in c(1:dim(ME[[1]])[[2]])){
		barplot(ME[[1]][o,i],col=col_geno[!exc][o],main=n[i])
		barplot(ME[[1]][o,i],density=d[!exc][o],col="white",add=TRUE)
	}
	dev.off()
}

m_names <- dimnames(ME[[1]])[[2]]
m <- apply(kME,1,which.max)
net <- cbind(dimnames(kME)[[1]],ens[dimnames(kME)[[1]],6],m_names[m],kME[cbind(1:dim(kME)[[1]],m)])

##Heatmaps early v late
hclust_a <- function(x){return(hclust(x,method="ward.D"))}
a <- cor(exp_matrix[cv>.05,b_dur_num[!exc]==0])
heatmap.2(a,trace="none",hclustfun=hclust_a,RowSideColors=col_geno[!exc][b_dur_num[!exc]==0],ColSideColors=col_tx[!exc][b_dur_num[!exc]==0])
a <- cor(exp_matrix[cv>.05,b_dur_num[!exc]==1])
heatmap.2(a,trace="none",hclustfun=hclust_a,RowSideColors=col_geno[!exc][b_dur_num[!exc]==1],ColSideColors=col_tx[!exc][b_dur_num[!exc]==1])

##LIMMA
exp_dat_early <- exp_matrix[,x]
y <- paste(b_geno[!exc][x],b_tx[!exc][x],sep="_")
f <- factor(y)
early_design <- model.matrix(~0+f)
early_fit <- lmFit(exp_dat_early,early_design)
cont_mat <- makeContrasts(vc=fNULL_VC-fWT_VC,rapa=fNULL_RAPA-fWT_RAPA,levels=early_design)
early_fit2 <- contrasts.fit(early_fit,cont_mat)
early_fit2 <- eBayes(early_fit2)
early_vc <- topTable(early_fit2,coef=1,n=Inf)
early_rapa <- topTable(early_fit2,coef=2,n=Inf)

exp_dat_late <- exp_matrix[,!x]
y <- paste(b_geno[!exc][!x],b_tx[!exc][!x],sep="_")
f <- factor(y)
late_design <- model.matrix(~0+f)
late_fit <- lmFit(exp_dat_late,late_design)
cont_mat <- makeContrasts(vc=fNULL_VC-fWT_VC,rapa=fNULL_RAPA-fWT_RAPA,levels=late_design)
late_fit2 <- contrasts.fit(late_fit,cont_mat)
late_fit2 <- eBayes(late_fit2)
late_vc <- topTable(late_fit2,coef=1,n=Inf)
late_rapa <- topTable(late_fit2,coef=2,n=Inf)

###WGCNA Null samples
x <- b_geno_num[!exc] == 2
y <- apply(exp_matrix[,x],1,sum) > 0
cv <- apply(exp_matrix[y,x],1,var)/apply(exp_matrix[y,x],1,mean)
pickSoftThreshold(t(net_dat[cv>0.05,]))
adj_mat <- cor(t(net_dat[cv>0.05,]))
a <- adj_mat < 0
adj_mat[a] <- 0
tom_mat <- TOMdist(adj_mat^6)

###
ME9_ens <- rownames(kME)[mods==9]
ME9_dat <- t(exp_matrix[ME9_ens,])
adjmat <- cor(ME9_dat)
a <- adjmat < 0
adjmat[a] <- 0
tommat <- TOMdist(adjmat^6)

cyto_gene_1 <- list()
cyto_gene_2 <- list()
cyto_weight <- list()

pb <- txtProgressBar(min=0,max=578,style=3)
for(i in c(1:578)){
	k <- c((i+1):579)
	cyto_gene_1[[i]] <- rep("",length(k))
	cyto_gene_2[[i]] <- rep("",length(k))
	cyto_weight[[i]] <- rep(NA,length(k))
	for(j in k){
		cyto_gene_1[[i]][j-i] <- ens[ME9_ens[i],6]
		cyto_gene_2[[i]][j-i] <- ens[ME9_ens[j],6]
		cyto_weight[[i]][j-i] <- tommat[i,j]
	}
	setTxtProgressBar(pb,i)
}